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Abstract. - Motivated by the dynamics of resonant neurons we consider a different iable, non- 
Markovian random process x(t) and particularly the time after which it will reach a certain level 
Xb- The probability density of this first passage time is expressed as infinite series of integrals 
over joint probability densities of x and its velocity x. Approximating higher order terms of 
this series through the lower order ones leads to closed expressions in the cases of vanishing and 
moderate correlations between subsequent crossings of xt- For a linear oscillator driven by white 
or coloured Gaussian noise, which models a resonant neuron, we show that these approximations 
reproduce the complex structures of the first passage time densities characteristic for the 
underdamped dynamics, where Markovian approximations (giving monotonous first passage 
time distribution) fail. 



Introduction. - The first passage time (FPT) of a stochastic process x(t) starting at t = 
from a given initial value within an a priory prescribed domain O of its state space is the 
time T when x(t) leaves fi for the first time. This concept was originally introduced by E. 
Schrodingcr when discussing behaviour of Brownian particles in external fields [1]. A large 
variety of problems ranging from noise in vacuum tubes, chemical reactions and nucleation [2] 
to stochastic resonance [3] , behaviour of neurons [4] , and risk management in finance [5] can 
be reduced to FPT problems. 

We will assign T(T) to be the flux through the boundary of f2 at time T, i.e. the probability 
density of the first passage time. Approaches to find T(T) are typically based either on the 
Fokker-Planck equation with an absorbing boundary [6] or on the renewal theory [7] . Despite 
the long history, explicit expressions for the FPT density are known only for a few cases. These 
include overdamped particles in the force free case, with time-independent constant forces and 
linear forces under influence of white noise [4, 8, 9, 10], as well as the case of a constant force 
under coloured noise [11]. Reasonable approximations exist for a few nonlinear forces [12, 13]. 

If the relaxation time r re i of the system to a quasistationary distribution in Q is much smaller 
than the typical first passage time, the escape from Q occurs from this quasiequilibrium state 
independent of the initial condition. Escapes occur with a constant rate inversely proportional 

Typeset using EURO-TeX 



2 



EUROPHYSICS LETTERS 



to the mean FPT. This is the case for many chemical reactions and nucleation processes [2] 
as well as for the leaky integrate-and-fire model of a neuron [14]. On times T exceeding r re i 
the FPT probability density T{T) decays exponentially. 

If the scale separation between relaxation and escape times does not hold, the dependence 
on the initial conditions gets crucial [20]. This situation is found in resonant neurons [15, 
16, 17, 18, 19] where x, here the voltage variable, exhibits damped subthreshold oscillations 
around the rest state being the attractor for deterministic dynamics of the system. If x(t) 
exceeds the excitation threshold Xb, a new spike begins. After spiking the voltage is reset to 
a fixed value x far from the rest state and can reach x h again prior to relaxation. T(T) 
giving the probability density function (PDF) of intervals between two consecutive spikes 
strongly deviates from an exponential. This is the situation we have in mind when developing 
the approach for obtaining T(T) for a non-Markovian differentiable random processes x{t) 
starting for t = at x(0) = x < Xb- 

Our approach is based on a series expansion for T(T) which is known for several decades 
[23]. But it was never used for explicit calculations. The approach is based on the theory 
of level crossings originally put forward by S.O. Rice [21]. A generalisation of his approach 
(based on what is called the Wiener- Rice series [22]) was used by Stratonovich to estimate the 
mean time spent by a stochastic process above the given level Xb [10]. The exact expression 
for T(T) is analogous to the Wiener-Rice series, but corresponds to the case when the initial 
value x(0) differs from the threshold x&. 

We first give an elementary derivation of this result and show then how this series can be used 
to obtain analytical expressions based on decoupling approximations. Explicit calculations are 
performed for an underdamped harmonic oscillator driven by white or coloured Gaussian noise, 
i.e. for a resonate-and-fire neuron with subthreshold oscillations and reset. The linearity of 
the model simplifies calculations but is not crucial for the approach. 

Level crossings and first passage times. - Let us first discuss the probability n\(xb, t\{xo}, 0)dt 
that a differentiable random process x(t) starting from a fully defined initial condition 
{xo} = {x{0), x(0), ...} (corresponding to its Markovian embedding) with x(0) < Xb crosses the 
level Xb between t and t+dt with positive velocity v(t) = x(t) > (i.e. performs an upcrossing). 
The upcrossing can only occur for such x(t) that Xb — vdt < x(t) < Xb- The probability 
of this is J^_ vdt P(x,v,t\{xo},0)dx = \v\P(xb, v, t\{x n }, 0)dt, where P(x, v, t\{x }, 0) is the 
joint PDF for x and v under given initial conditions. Integration over all v > then gives 
ni(xb, t\{xo}, 0) = J Q vP(xb, v, t\{xo}, 0)dv. The joint probabilities for multiple upcrossings 
n p (t p , . . . , ti)dt p . . . dt\ in each of p time intervals (ti, t\ + dt\), . . . (t p , t p + dt p ) follow in the 
same way: 

n p (t p ,...,ti)= dv p ... dv 1 v p ...viP(xb,v p ,t p ;...;xb,vi,ti\{x },0). (1) 
Jo Jo 

The initial conditions and Xb will be omitted in what follows. 

The function !F(T) is given by the fraction of all trajectories x(t) which perform the first 
upcrossing of Xb at time T. All such trajectories are accounted for in n\{T). However, n\{T) 
also considers trajectories which might have an earlier upcrossing at some < t\ < T. Since 
they should not contribute to T(T), their fraction should be subtracted from n\{T) by taking 
ni(T) — J Q n 2 (T,ti)dti. This excludes all trajectories which cross x b exactly twice until 
T. However the trajectories crossing Xb three times, i.e at T and at two earlier moments 
ti < T, (i = 1, 2), are not correctly accounted. Each such trajectory is added once in ni(T), but 
subtracted two times in J Q T n 2 (T, ti)dti, since this term accounts for the pairs of upcrossings 
at (T,ti) with i = 1,2. To account for this we have to add once the amount of trajectories 



T. VERECHTCHAGUINA et al.\ : FIRST PASSAGE TIME DENSITIES 



3 



with three upcrossings again: m (T) — J Q n 2 (T, ti)dti + ^ J Q J Q n 3 (T, t 2 , ti)dt 2 dti (the factor 
1/2! in the last term accounts for permutations of ti). 

Generally, if a trajectory crosses Xb at time T and at N earlier times ti<T,i = l...N, 
then in ^ J ... J n p+ i(T, t p , . . . , t\)dt p . . . dti it is accounted for exactly Cj^ times (Cj^ stands 

for the number of combinations). Note, that X^=o( — -0 P ^iv = (1 — 1) W = Thus in the 
alternating sum of this kind containing N + 1 terms, all trajectories crossing Xb at time T 
and having 1, 2, . . . N additional upcrossings are excluded. Extending the sum to infinity we 
exclude all superfluous trajectories. Only trajectories remain for which the upcrossing at time 
T was the first one. Thus the expression for the first passage time probability density reads: 

°° (-l) p f T f T 

n p+ i(T,t p , . . . (2) 

p=o p - Jo Jo 

In Rcf. [7, 10] it was shown that !F(T) can be equivalently expressed in terms of the 
correlation functions between upcrossings (cross-cumulants) g p (t p , . . . , t\). They are related 
to the joint densities n p (t p , . . . , t\) via 

gi(ti) = rai(ii), 

92(t2,h) = n 2 (t 2 ,t 1 ) - ni(ii)m(i 2 ), (3) 
ff3(*3,*2,*i) = n3(<3,*2,ti) - 3{ni(ti)n 2 (i3,i2)} ;s + 2ni(ii)ni(t 2 )ni(t 3 ), 

Here {. . .} s denotes the symmetrisation of the expression in the brackets with respect to all 
permutations of its arguments. The expression for T{T) in terms of correlation functions reads 
[10]: 

T(T) = S'(T)e- s ^ (4) 

with 

°° (— l) p f T f T 
S(T) = -J2 y —^ dt p ... dt l9p {t p ,...,h)- (5) 
p =i p - Jo Jo 

The function S'(T) can be interpreted as a the time-dependent escape rate. 

We emphasise that expressions eqs.(2),(4) and (5) do hold for all differentiable random 
processes (ones whose velocity v(t) is defined at any time). 

Decoupling approximations for the FPT density. Dealing with infinite series useful 
approximations can either be based on the truncation of the series after several first terms 
calculated exactly, or on expressing the higher order terms approximately through the lower 
order ones what might lead to a closed analytical form. Truncation approximations for eq.(2) 
are not normalised and hold on short time scales only diverging at longer times (due to the 
miscount of trajectories with several upcrossings). We discuss here the approximations of 
the second type which are normalised and can be used in the whole time domain. Each 
such approximation is based on a subsummation in eq.(5) for S(T). Note, only expressions 
guaranteeing positive rates S(T) are allowed. 

The simplest approximation is based on neglecting all terms in eq.(5) except for the first 
one. It leads to 



S(T) = [ m(t)dt, 

Jo 



(6) 



4 



EUROPHYSICS LETTERS 



where n\{t) is given by eq.(l) with p = 1. This corresponds to neglecting all correlations 
between upcrossings of the level Xb, i.e. to factorisation of n p+ i(T, t p , . . . , ti) into a product 
of one-point densities ni(T)n\(t p ) . . . n\(ti) in cq.(2). Then the series, eq.(2) sums up into 

J-(T) Ri ni(T) exp ^— ni(t)dtj , which is equivalent to cq.(6). The approximation will be 

refereed to as a Hertz approximation since the form of F(T) resembles the Hertz distribution 
[24]. 

The second order approximation expresses all higher order correlation functions through 
the first and the second ones. It was used by Stratonovich [10] in the context of peak duration. 
The first and the second correlation functions are taken exactly, and the higher ones are 
approximated by the combinations of these two. For p > 2 one thus has 

g p (t p , ...,*!)« {-\y-\p - l)!ni(t p ) . . . ni(ti){i2(ti, t 2 )R(h, t 3 ) . . . R(h,t p )} s (7) 

with the correlation coefficient 

n\{ti)n\{tj) 

The correlation coefficient i?(ii,i 2 ) is equal to unity for t\ = ti and vanishes for large values 
of |ti — ti\. Substitution of eq.(7) into eq.(5) delivers then the Stratonovich approximation for 
T{T) with the time-dependent escape rate being 



= -/ »i(t) 
Jo 



In 



l - f T i?(i, t')m{t')dt' 



S{T) = - m{t) — A -dt. (9) 

/„ R(t,t') ni (t')dt' 

Noise driven harmonic oscillator. - Let us now illustrate the method by applying it to the 
coordinate x(t) of a harmonic oscillator driven by a Gaussian noise n{t) 

x = v, v = — jv — lUqX + rj(t). (10) 

Two cases will be considered: (i) white noise, r)(t) = \/2D^(t), and (ii) Ornstein-Uhlenbeck 
noise, f\ = —T~ x r\ + \> / 2DT~ 1 £ i (t), with ^(t) being white Gaussian noise of intensity 1. Eq.(10) 
with boundary at Xb and reset at x(0) = xq,x(0) = vq is relevant for a modelling of neuronal 
dynamics. In the undcrdamped regime 7 < 2ujq it describes an excitable dynamics with 
damped subthreshold oscillations characteristic for resonant neurons [15, 16, 17, 18, 19]. In 
the overdamped regime 7 > 2luq it describes the behaviour of nonresonant neurons [15, 25]. 
In our calculations we fix wo = 1, take initial conditions to be xq = — 1,vq = 0, and set the 
absorbing boundary at the threshold Xb = 1. 

Our basic model has an advantage that all transition probability densities in eq.(l) are 
Gaussian and can be expressed through the correlation functions of their arguments, which 
are easily obtained using the spectral representation [6] . Then n\ (T) is obtained in closed ana- 
lytical form; the joint densities of multiple upcrossings are readily obtained through numerical 
evaluation of integrals in eq.(l). 

The Hertz and the Stratonovich approximations cqs.(6) and (9) hold if all correlations decay 
considerably within the typical time interval between upcrossings. The decay of correlations 
is described by the relaxation time T re i = 2/7 of the process. The mean interval between 
two successive upcrossings Tr is the inverse stationary frequency of upcrossings 1/Tr = uq = 
lim^oo n\{t) [21]. n is known as the Rice frequency and is given by 
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Fig. 1. - FPT probability density for harmonic oscillator driven by Gaussian white noise. Simulation 
results are shown with a gray line, Hertz approximation with a black dashed line, and Stratonovich 
approximation with a black solid line. Note the logarithmic scale in T. The insets show the same 
curves on the logarithmic scale in T(T). The parameters are: ujo = l,a;o = — l,«o = 0, xt, = 1, 
a) 7 = 0.8, D = 0.1, TVei = 2.5, T R = 343, b) 7 = 0.8, D = 0.44, r rel = 2.5, T R = 15.6, c) 
7 = 0.08, D = 0.01, r rel = 25, T R = 343, d) 7 = 10.0, D = 5.5, W0/7 = 0.1 



where r xx (f) = (x(t)x(O)) is the correlation function of the process. In our case Tr = 
27T(jj e~ /XbU ' o/ . Therefore, the Hertz approximation holds for T re \ -c 7^; for the Stratonovich 
approximation this condition can be weakened to r re ; < Tr arising from the condition that 
the argument of the logarithm in eq.(9) is positive. 

Let us first concentrate on the case (i) with white noise. In fig. 1 the FPT probability den- 
sity obtained from simulations is compared with the Hertz and Stratonovich approximations 
(eq.(6) and eq.(9)). The probability to reach xt is higher in the maxima of the subthreshold 
oscillations. The initial phase of these oscillations is fixed by initial conditions. Thus on 
shorter time scales T(T) shows the multiple peaks following with the frequency of damped 
oscillations ^ u>^ — 7 2 /4. On long times T 3> r re i the quasicquilibrium establishes and FPT 
densities decay exponentially (see insets in fig. 1). The number of visible peaks depends on 
the relation between r re i and the period of oscillations and is given by the number of periods 
elapsing before the quasiequilibrium is achieved. 

In fig. 1(a) the parameters are chosen to be 7 = 0.8, D = 0.1, corresponding to moderate 
friction and moderate noise intensity. For given parameter values r re i = 2.5 and Tr = 343, 
so that r re i <C Tr, both Hertz and Stratonovich approximations hold and reproduce well the 
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FPT density in the whole time domain. 

In the case of moderate friction and stronger noise the upcrossings become more frequent 
and Tr decreases. The FPT changes its form to practically monomodal with the only maximum 
and a small shoulder separating it from the tail. An example is given in fig. 1(b) with 
7 = 0.8, D = 0.44 which correspond to r re i = 2.5 and Tr = 15.6. The Stratonovich 
approximation complies very well with the simulations, while the Hertz approximation fails to 
reproduce the details of the distribution: It underestimates T(T) on short times, and shows 
slower exponential decay in the tail than the one observed in simulations (see the inset). 

Finally, for small friction and low noise the upcrossings are rare, but the relaxation time is 
large. The FPT probability density exhibits multiple decaying peaks. In fig. 1(c) 7 = 0.08, D = 
0.01 corresponding to r re ; = 25, Tr = 343. Again, the Stratonovich approximation performs 
well, while the Hertz approximation underestimates the first peak, overestimates all further 
peaks and decays in the tail faster than the simulated FPT density. 

Fig. 1(d) corresponds to the overdamped regime (7 > 2cjo) where the parameters are 
chosen to be 7 = 10.0,1? = 5.5 corresponding to ojq/j = 0.1. The condition r re i < Tr 
is always fulfilled in the overdamped case. However, with increasing friction the process 
x(t) approaches the Markovian (diffusion) one, for which the pattern of upcrossings is very 
inhomogeneous. The successive crossings are strongly clustered [10]. The upcrossings within 
a cluster are not independent even if their mean density Hq is low. This limits the accuracy of 
our approximations: The Stratonovich approximation starts to be inaccurate, and the Hertz 
approximation fails. 

For T large T(T) decays exponentially, T(T) oc exp(— kT). The decrement of this decay 
is obtained from long time asymptotics: kT = limr^ 00 5(T). Thus, in the Hertz approx- 
imation eq.(6) one gets k = (1/T) limr^oo ni(t)dt = n T/T = n . The behaviour in 
the Stratonovich approximation eq.(9) is determined by limt,t'->oo R(t, t')m(t')dt' w noT cor 
with T cor given by T cor — lim f n °° R(t,t')dt' . Note that r cor is not necessary positive. In- 

serting this expression into eq.(9) and expanding the logarithm up to the second term we get 
k = no(l + \noT cor ) providing the second order correction to the previous expression. The 
value of r cor for the parameter set as in fig. 1(a) is r cor = —1 • 10~ 3 , for parameters as in fig. 
1(b) T cor = 5.25, and for parameters as in fig. 1(c) r cor = —396.7. The long time asymptotics 
obtained with these r cor values reproduce fairly well the decay patterns found numerically. 

Finally we consider our model (ii) having a higher dimension of its state space. In fig. 2 
we show the simulated FPT probability density and the Hertz approximation for the system 
eq.(10) driven by the Ornstcin-Uhlenbcck noise for two different values of the correlation time 
t = 0.5 and r = 1.0. In both cases the Hertz approximation is absolutely sufficient. 

Conclusions. - For the case of moderate friction and moderate noise intensity the Hertz 
approximation is absolutely sufficient. The Stratonovich approximation performs evenly well 
and does not lose accuracy for high noise intensity. The validity region of approximations 
covers all different types of subthreshold dynamics, and reproduces all qualitative different 
structures of the FPT PDF: from monomodal through bimodal to multimodal densities with 
decaying peaks. The approximations work for the systems of whatever dimension and are 
especially effective for the processes with narrow spectral density, exactly when Markovian 
approximations fail. 

We acknowledge financial support from the DFG by Graduierten-Kolleg 268 and the Bern- 
stein Center for Computational Neuroscience, Berlin and Sfb555. 
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Fig. 2. - Simulation results and the Hertz approximation for the FPT probability density for harmonic 
oscillator driven by Ornstein-Uhlenbeck noise with correlation times r = 0.5 (upper curves) and 
r = 1.0 (lower curves) and other parameters as in fig. 1(a). Note the logarithmic scale in T. The 
inset shows the same curves on the logarithmic scale in T(T) . 
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